function y_=simult_amd(M_,options_,y0,dr,ex_,iorder)%,aux)
% Simulates the model using a perturbation approach, given the path for the exogenous variables and the
% decision rules.
%
% INPUTS
%    M_       [struct]   model
%    options_ [struct]   options
%    y0       [double]   n*1 vector, initial value (n is the number of declared endogenous variables plus the number
%                        of auxilliary variables for lags and leads); must be in declaration order, i.e. as in M_.endo_names
%    dr       [struct]   matlab's structure where the reduced form solution of the model is stored.
%    ex_      [double]   T*q matrix of innovations.
%    iorder   [integer]  order of the taylor approximation.
%
% OUTPUTS
%    y_       [double]   n*(T+1) time series for the endogenous variables.
%
% SPECIAL REQUIREMENTS
%    none

% Copyright (C) 2001-2019 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.



oo_.dr = dr;

iter = size(ex_,1);
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;

y_ = zeros(size(y0,1),iter+M_.maximum_lag);
y_(:,1) = y0;

if options_.loglinear && ~options_.logged_steady_state
    dr.ys=log(dr.ys);
end

if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if k_order_pert is not used or if we do not use Dynare++ with k_order_pert
    if iorder==1
        y_(:,1) = y_(:,1)-dr.ys;
    end
end

if options_.k_order_solver && ~options_.pruning % Call dynare++ routines.
    ex_ = [zeros(M_.maximum_lag,M_.exo_nbr); ex_];
    [err, y_] = dynare_simul_(options_.order,M_.nstatic,M_.npred,M_.nboth,M_.nfwrd,exo_nbr, ...
        y_(dr.order_var,1),ex_',M_.Sigma_e,options_.DynareRandomStreams.seed, ...
        dr.ys(dr.order_var),dr);
    mexErrCheck('dynare_simul_', err);
    y_(dr.order_var,:) = y_;
else
    if options_.block
        if M_.maximum_lag > 0
            k2 = dr.state_var;
        else
            k2 = [];
        end
        order_var = 1:endo_nbr;
        dr.order_var = order_var;
    else
        k2 = dr.kstate(find(dr.kstate(:,2) <= M_.maximum_lag+1),[1 2]);
        k2 = k2(:,1)+(M_.maximum_lag+1-k2(:,2))*endo_nbr;
        order_var = dr.order_var;
    end
    
    switch iorder
        case 1
            if isempty(dr.ghu)% For (linearized) deterministic models.
                for i = 2:iter+M_.maximum_lag
                    yhat = y_(order_var(k2),i-1);
                    y_(order_var,i) = dr.ghx*yhat;
                end
            elseif isempty(dr.ghx)% For (linearized) purely forward variables (no state variables).
                y_(dr.order_var,:) = dr.ghu*transpose(ex_);
            else
                epsilon = dr.ghu*transpose(ex_);
                for i = 2:iter+M_.maximum_lag
                    yhat = y_(order_var(k2),i-1);
                    y_(order_var,i) = dr.ghx*yhat + epsilon(:,i-1);
                end
            end
            y_ = bsxfun(@plus,y_,dr.ys);
        case 2
            constant = dr.ys(order_var)+.5*dr.ghs2;
            if options_.pruning
                y__ = y0;
                for i = 2:iter+M_.maximum_lag
                    yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
                    yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
                    epsilon = ex_(i-1,:)';
                    [abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1);
                     
                    [abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
                     
                    [abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon);
                     
                    y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
                        + abcOut1 + abcOut2 + abcOut3;
                    y__(order_var) = dr.ys(order_var) + dr.ghx*yhat1 + dr.ghu*epsilon;
                end
            else
                for i = 2:iter+M_.maximum_lag
                    yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
                    epsilon = ex_(i-1,:)';
                    [abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat);
                     
                    [abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon);
                     
                    [abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon);
                     
                    y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
                        + abcOut1 + abcOut2 + abcOut3;
                end
            end
        case 3
            % only with pruning
            % the third moments of the shocks are assumed null. We don't have
            % an interface for specifying them
            ghx = dr.ghx;
            ghu = dr.ghu;
            ghxx = dr.ghxx;
            ghxu = dr.ghxu;
            ghuu = dr.ghuu;
            ghs2 = dr.ghs2;
            ghxxx = dr.ghxxx;
            ghxxu = dr.ghxxu;
            ghxuu = dr.ghxuu;
            ghuuu = dr.ghuuu;
            ghxss = dr.ghxss;
            ghuss = dr.ghuss;
            

            
            nspred = M_.nspred;
            ipred = M_.nstatic+(1:nspred);
            %construction follows Andreasen et al (2013), Technical
            %Appendix, Formulas (65) and (66)
            %split into first, second, and third order terms
            yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
            yhat2 = zeros(nspred,1);
            yhat3 = zeros(nspred,1);
            for i=2:iter+M_.maximum_lag
                u = ex_(i-1,:)';
                %construct terms of order 2 from second order part, based
                %on linear component yhat1
                [gyy] = A_times_B_kronecker_C(ghxx,yhat1);
                % 
                [guu] = A_times_B_kronecker_C(ghuu,u);
                % 
                [gyu] = A_times_B_kronecker_C(ghxu,yhat1,u);
                % 
                %construct terms of order 3 from second order part, based
                %on order 2 component yhat2
                [gyy12] = A_times_B_kronecker_C(ghxx,yhat1,yhat2);
                 
                [gy2u] = A_times_B_kronecker_C(ghxu,yhat2,u);
                 
                %construct terms of order 3, all based on first order component yhat1
                y2a = kron(yhat1,yhat1);
                [gyyy] = A_times_B_kronecker_C(ghxxx,y2a,yhat1);
                 
                u2a = kron(u,u);
                [guuu] = A_times_B_kronecker_C(ghuuu,u2a,u);
                 
                yu = kron(yhat1,u);
                [gyyu] = A_times_B_kronecker_C(ghxxu,yhat1,yu);
                 
                [gyuu] = A_times_B_kronecker_C(ghxuu,yu,u);
                 
                %add all terms of order 3, linear component based on third
                %order yhat3
                yhat3 = ghx*yhat3 +gyy12 ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
                    + gy2u ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
                    + 1/6*(gyyy + guuu + 3*(gyyu + gyuu +  ghxss*yhat1 + ghuss*u)); %note: s is treated as variable, thus xss and uss are third order
                yhat2 = ghx*yhat2 + 1/2*(gyy + guu + 2*gyu + ghs2);
                yhat1 = ghx*yhat1 + ghu*u;
                y_(order_var,i) = dr.ys(order_var)+yhat1 + yhat2 + yhat3; %combine terms again
                yhat1 = yhat1(ipred);
                yhat2 = yhat2(ipred);
                yhat3 = yhat3(ipred);
            end
        case 4
            % only with pruning
            % the third moments of the shocks are assumed null. We don't have
            % an interface for specifying them
            ghx = dr.ghx;
            ghu = dr.ghu;
            ghxx = dr.ghxx;
            ghxu = dr.ghxu;
            ghuu = dr.ghuu;
            ghs2 = dr.ghs2;
            ghxxx = dr.ghxxx;
            ghxxu = dr.ghxxu;
            ghxuu = dr.ghxuu;
            ghuuu = dr.ghuuu;
            ghxss = dr.ghxss;
            ghuss = dr.ghuss;
            
            %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Fourth order AMD 1/10/2021
            
            G4 = oo_.dr.g_4;  % This is the folded g_4 matrix w
            G3 = oo_.dr.g_3;
            G2 = oo_.dr.g_2;
            G1 = oo_.dr.g_1;
            n=size(G1,2);
            
            G2_new = unfoldg2(G2,n);
            
            G3_new = unfoldg3(G3,n);

            nn = size(oo_.dr.g_1,2); % This gets number of state vars + shocks
            n_x = M_.npred+M_.nboth;
            n_st = size(oo_.dr.g_1,2);
            n_exo = n_st - n_x;
            G4_new = unfoldg4(G4,nn); % This calls a function that unfolds g_4
            
            %Fourth order components that show up in the second order
            G4_ghxx=[];
            for i=1:n_x
                G4_ghxx = [G4_ghxx G2_new( :, 1+n_st*(i-1) : 1+n_st*(i-1)+n_x-1)];
            end
            ghs2xx = (G4_ghxx - oo_.dr.ghxx*0.5)/(6/24);
            %At the fourth order approximation  G_2 = (1/2)*ghxx + (6/24)*ghs2xx
            
            %Isolating the sigma_a x sigma_a column
            
%             for i=1:49
%                 if i==aux   
%                 else
%                     ghs2xx(:,i) = zeros(47,1);
%                 end
%             end
            
            G4_ghxu=[];
            for i=1:n_x
                G4_ghxu = [G4_ghxu G2_new( :, 1+n_x+n_st*(i-1) : 1+n_x+n_st*(i-1)+n_exo-1)];
            end
            %At the third order approximation  2*G_2 = (1/2)*2*ghxu
            %Because I can show that G_2 = 1/2*ghxu.
            %At the fourth order approximation 2*G_2 = (1/2)*2*ghxu + 12*(1/24)*ghs2xu
            ghs2xu = (2*G4_ghxu - 2*oo_.dr.ghxu*0.5)/(12/24);
            
                    %Isolating the sigma_a x sigma_a column
%             for i=1:28
%                 if i==aux   
%                 else
%                     ghs2xu(:,i) = zeros(47,1);
%                 end
%             end    
            
            
            %oo_.dr.g_0 = 0.5*dr.ghs2 + (1/24)*ghs4
            %The ghs3 terms are zero
            ghs4 = (oo_.dr.g_0 - 0.5*dr.ghs2)/(1/24);
            
            
            G4_ghuu=[];
            for i=1:n_exo
                G4_ghuu = [G4_ghuu G2_new( :, n_x*n_st + n_x + 1 +n_st*(i-1) :n_x*n_st + n_x + 1 +n_exo-1 +n_st*(i-1))];
            end
            %At the fourth order approximation  G_2 = (1/2)*ghuu + (6/24)*ghs2uu
            ghs2uu = (G4_ghuu - oo_.dr.ghuu*0.5)/(6/24);
            
            
            
            G4_ghxxxx=[];
            for iii =  1:n_x
                for ii = 1:n_x
                    for i=1:n_x
                        G4_ghxxxx = [G4_ghxxxx G4_new( :, 1	+									n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	:	1	+									n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	+	n_x -1)];
                    end
                end
            end
            % G_4  = (1/24)*oo_.dr.ghxxxx
            ghxxxx = 24*G4_ghxxxx;
            
            G4_ghxxxu=[];
            for iii =  1:n_x
                for ii = 1:n_x
                    for i=1:n_x
                        G4_ghxxxu = [G4_ghxxxu G4_new( :, 1	+	n_x	+		+		+		+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	:	1	+	n_x	+		+		+		+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	+	n_exo-1)];
                    end
                end
            end
            ghxxxu = 24*G4_ghxxxu;
            
            G4_ghxxuu=[];
            for iii =  1:n_x
                for ii = 1:n_x
                    for i=1:n_exo
                        G4_ghxxuu = [G4_ghxxuu G4_new( :, 1	+	n_x	+	n_x*n_st	+		+		+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	:	1	+	n_x	+	n_x*n_st	+		+		+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	+	n_exo-1)];
                    end
                end
            end
            ghxxuu = 24*G4_ghxxuu;
            
            G4_ghxuuu=[];
%           disp(num2str(aux));
            for iii =  1:n_x
                for ii = 1:n_exo
                    for i=1:n_exo
                        G4_ghxuuu = [G4_ghxuuu G4_new( :, 1	+	n_x	+	n_x*n_st	+	n_x*n_st^2	+		+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	:	1	+	n_x	+	n_x*n_st	+	n_x*n_st^2	+		+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	+	n_exo-1)];
                    end
                end
            end     
            ghxuuu = 24*G4_ghxuuu;
            
            G4_ghuuuu=[];
            for iii =  1:n_exo
                for ii = 1:n_exo
                    for i=1:n_exo
                        G4_ghuuuu = [G4_ghuuuu G4_new( :, 1	+	n_x	+	n_x*n_st	+	n_x*n_st^2	+	n_x*n_st^3	+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	:	1	+	n_x	+	n_x*n_st	+	n_x*n_st^2	+	n_x*n_st^3	+	n_st*(i-1)	+	n_st^2*(ii-1)	+	n_st^3*(iii-1)	+	n_exo-1)];
                    end
                end
            end
            ghuuuu = 24*G4_ghuuuu;
            %%
            
            threads = 1;% options_.threads.kronecker.A_times_B_kronecker_C;
            nspred = M_.nspred;
            ipred = M_.nstatic+(1:nspred);
            %construction follows Andreasen et al (2013), Technical
            %Appendix, Formulas (65) and (66)
            %split into first, second, and third order terms
            yhat1 = y0(order_var(k2))-dr.ys(order_var(k2));
            yhat2 = zeros(nspred,1);
            yhat3 = zeros(nspred,1);
            yhat4 = zeros(nspred,1);
            for i=2:iter+M_.maximum_lag
                u = ex_(i-1,:)';
                %construct terms of order 2 from second order part, based
                %on linear component yhat1
                [gyy] = A_times_B_kronecker_C(ghxx,yhat1);
                 
                [guu] = A_times_B_kronecker_C(ghuu,u);
                 
                [gyu] = A_times_B_kronecker_C(ghxu,yhat1,u);
                 
                
                
                %construct terms of order 3 from second order part, based
                %on order 2 component yhat2
                
                [gyy12] = A_times_B_kronecker_C(ghxx,yhat1,yhat2);
                 
                [gy2u] = A_times_B_kronecker_C(ghxu,yhat2,u);
                 
                
                
                %construct terms of order 3, all based on first order component yhat1
                
                y2a = kron(yhat1,yhat1);
                [gyyy] = A_times_B_kronecker_C(ghxxx,y2a,yhat1 );
                 
                u2a = kron(u,u);
                [guuu] = A_times_B_kronecker_C(ghuuu,u2a,u );
                 
                yu = kron(yhat1,u);
                [gyyu] = A_times_B_kronecker_C(ghxxu,yhat1,yu );
                 
                [gyuu] = A_times_B_kronecker_C(ghxuu,yu,u );
                 
                
                
                %construct terms of order 4, all already
                
                
                %add all terms of order 4, linear component based on fourth
                %order yhat4
                
                [g_hxx_yhat1_yhat3] = A_times_B_kronecker_C(ghxx,yhat1,yhat3 );
                 
                [g_hxx_yhat2_yhat2] = A_times_B_kronecker_C(ghxx,yhat2,yhat2 );
                 
                [g_hxu_yhat3_u] = A_times_B_kronecker_C(ghxu,yhat3,u );
                 
                
                y1y1 = kron(yhat1,yhat1);
                
                [g_hxxx_yhat1_yhat1_yhat2] = A_times_B_kronecker_C(ghxxx,y1y1,yhat2 );
                 
                
                y1y2 = kron(yhat1,yhat2);
                [g_hxxu_yhat1_yhat2_u] = A_times_B_kronecker_C(ghxxu,y1y2,u );
                 
                
                y2u = kron(yhat2,u);
                [g_hxuu_yhat2_u_u] = A_times_B_kronecker_C(ghxuu,y2u,u);
                 
                
                y1y1 = kron(yhat1,yhat1);
                y1y1y1 = kron(y1y1,yhat1);
                [g_hxxxx_yhat1_yhat1_yhat1_yhat1] = A_times_B_kronecker_C(ghxxxx,y1y1y1,yhat1 );
                 
                
                [g_hxxxu_yhat1_yhat1_yhat1_u] = A_times_B_kronecker_C(ghxxxu,y1y1y1,u );
                 
                
                y1y1u = kron(y1y1,u);
                [g_hxxuu_yhat1_yhat1_u_u] = A_times_B_kronecker_C(ghxxuu,y1y1u,u );
                 
                              
                y1u = kron(yhat1,u);
                y1uu = kron(y1u,u);
                [g_hxuuu_yhat1_u_u_u] = A_times_B_kronecker_C(ghxuuu,y1uu,u );
                 
                
                uu = kron(u,u);
                uuu = kron(uu,u);
                [g_huuuu_u_u_u_u] = A_times_B_kronecker_C(ghuuuu,uuu,u );
                 
                
                [g_hs2xx_yhat1_yhat1] = A_times_B_kronecker_C(ghs2xx,yhat1,yhat1 );
                 
                
                [g_hs2uu_u_u] = A_times_B_kronecker_C(ghs2uu,u,u );
                 
                
                [g_hs2xu_yhat1_u] = A_times_B_kronecker_C(ghs2xu,yhat1,u );
                 

                
                yhat4 = ghx*yhat4 + ...
                        (1/2) * ( 2* g_hxx_yhat1_yhat3 + g_hxx_yhat2_yhat2 + 2*g_hxu_yhat3_u) + ...
                        (1/6) * ( 3* g_hxxx_yhat1_yhat1_yhat2 + ...
                                  3* g_hxxu_yhat1_yhat2_u + ...
                                  3* g_hxuu_yhat2_u_u + ...
                                  3* ghxss*yhat2) + ...
                        (1/24)* ( 1*g_hxxxx_yhat1_yhat1_yhat1_yhat1 + ...
                                  1*4* g_hxxxu_yhat1_yhat1_yhat1_u + ...
                                  1*6* g_hxxuu_yhat1_yhat1_u_u + ...
                                  1*4* g_hxuuu_yhat1_u_u_u + ...
                                  1*g_huuuu_u_u_u_u + ...
                                  1*6*g_hs2xx_yhat1_yhat1 + ...
                                  1*6*g_hs2uu_u_u + ...
                                  1*12*g_hs2xu_yhat1_u + ...
                                  1*ghs4);
                % ghx*yhat4 +
                %
                % (1/2)*( g_hxx * (2*kron(yhat1,yhat3) + kron(yhat2,yhat2)) +
                % 2*g_hxu*kron(yhat3,u))
                % +
                % (1/6)*( g_hxxx*(3*kron(yhat1,yhat1,yhat2)) +
                % 3*g_hxxu*(2*kron(yhat1,yhat2,u) +
                % 3*g_hxuu*kron(yhat2,u,u) +
                % 3*ghxss*yhat2))
                % +
                % (1/24)( g_hxxxx*kron(yhat1,yhat1,yhat1,yhat1) +
                % 4*g_hxxxu*(yhat1,yhat1,yhat1,u) +
                % 6*g_hxxuu*kron(yhat1,yhat1,u,u) +
                % 4*g_hxuuu*kron(yhat1,u,u,u) + g_huuuu*kron(u,u,u,u) +
                %
                % 6*g_hs2xx*kron(yhat1,yhat1) +
                % g_hs2xx =  g_hxx solved at fourth order minus g_hxx at second order
                % ghs2xx = (ghxx - oo_.dr.ghxx*0.5)/(6/24);
                %
                % 6*g_hs2uu*kron(u,u) +
                % g_huu = g_huu solved at fourth order minus g_
                % ghs2uu = (ghuu - oo_.dr.ghxuu*0.5)/(6/24);
                %
                % 12*g_hs2xu*kron(yhat1,u) +
                % ghs2xu = (2*ghxu - 2*oo_.dr.ghxu*0.5)/(12/24);
                %
                % 4*g_hs3u*u +
                % This is practically zero
                % 4*g_hs3x*yhat1 +
                % This is practically zero
                % g_hs4 )
                %
                % ghs4 = (oo_.dr.g_0 - 0.5*dr.ghs2)/(1/24);

                %add all terms of order 3, linear component based on third
                %order yhat3
                
%                 if i==5000
%                     i
%                 end
                
                yhat3 = ghx*yhat3 +gyy12 ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
                    + gy2u ... % prefactor is 1/2*2=1, see (65) Appendix Andreasen et al.
                    + 1/6*(gyyy + guuu + 3*(gyyu + gyuu +  ghxss*yhat1 + ghuss*u)); %note: s is treated as variable, thus xss and uss are third order
                yhat2 = ghx*yhat2 + 1/2*(gyy + guu + 2*gyu + ghs2);
                yhat1 = ghx*yhat1 + ghu*u;
                y_(order_var,i) = dr.ys(order_var)+yhat1 + yhat2 + yhat3 + yhat4; %combine terms again
                yhat1 = yhat1(ipred);  %this converts yhats back to state variables
                yhat2 = yhat2(ipred);
                yhat3 = yhat3(ipred);
                yhat4 = yhat4(ipred);
                
            end
            disp('4th Order Pruning Complete');
        otherwise
            error(['pruning not available for order = ' int2str(options_.order)])
    end
end
